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^ Abstract 

<! 

_j The semi-analytical wall boundary conditions present a mathematically rigorous framework 

1-H to prescribe the influence of solid walls in SPH for fluid flows. In this paper they are inves- 

I I tigated with respect to the skew-adjoint property which implies exact energy conservation. 

C It will be shown that this property holds only in the limit of the continuous SPH approx- 

,^ imation, whereas in the discrete SPH formulation it is only approximately true, leading to 

i± numerical noise. This noise, interpreted as form of "turbulence", is treated using an addi- 

r—j tional volume diffusion term in the continuity equation which we show is equivalent to an 

j/5 approximate Riemann solver. Subsequently two extensions to the boundary conditions are 

.^ presented. The first dealing with a variable driving force when imposing a volume flux in 

^ a periodic flow and the second showing a generalization of the wall boundary condition to 

^ Robin type and arbitrary-order interpolation. Two modifications for free-surface flows are 

, ^ presented for the volume diffusion term as well as the wall boundary condition. In order to 

validate the theoretical constructs numerical experiments are performed showing that the 

►^ present volume flux term yields results with an error 5 orders of magnitude smaller then 

(N previous methods while the Robin boundary conditions are imposed correctly with an error 

^ depending on the order of the approximation. Furthermore, the proposed modifications for 

(T^ free-surface flows improve the behaviour at the intersection of free surface and wall as well as 

■^ prevent free-surface detachment when using the volume diffusion term. Finally, this paper 

O is concluded by a simulation of a dam break over a wedge demonstrating the improvements 
proposed in this paper. 
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1. Introduction 

Smoothed particle liydrodynamics (SPH) lias become one of the predominant meshless 
methods to solve the Navier-Stokes equations. An introduction to this method can be found 
in the reference paper by Monaghan [T| or the book by Violeau [2]. One of the greatest 
difficulties associated with this method is the imposition of boundary conditions. The semi- 
analytical wall boundary conditions introduced by Ferrand et al. [3J provide a consistent 
framework for imposing wall boundary conditions in both laminar and turbulent flows. Con- 
trary to previous approaches such as the Leonnard- Jones potential forces (see Monaghan |3]) 
or more advanced repulsive force models (see Monaghan and Kajtar [3]) they do not rely on 
additional modelling but instead are introduced naturally from the numerical treatment of 
the governing equations. Furthermore, they can deal with arbitrarily complex geometries in 
contrast to the approach using ghost particles (see e.g. Colagrossi and Landrini [6]). 
Despite its ability to predict flows with good accuracy, SPH with the semi-analytical wall 
boundary conditions still suffers from numerical drawbacks. This can be observed, for exam- 
ple, when considering the energy budget from an inviscid flow between two infinite plates. 
As the energy is not constant the skew-adjoint principle is violated thus warranting an in- 
vestigation of this property with respect to the semi-analytical wall boundary conditions. In 
addition to a theoretical analysis a numerical simulation will show stability issues associated 
with this property. This leads to the question of whether it is possible to counter the insta- 
bilities by some numerical procedure. If a diffusion term is chosen, as in this paper, then 
its properties and influence on the fluid flow need to be investigated. The question of how 
to impose a volume flux or speciflc boundary conditions are also closely linked to modelling 
wall boundary conditions. In free-surface flows certain difficulties can occur which are not 
apparent in confined flows. To solve these, it is necessary to revisit certain investigations by 
separating hydrostatic and dynamic pressure differences. 

In the following section the wall boundary conditions will be introduced in detail and will 
then be investigated with respect to the skew-adjoint property in Section [3j This will be 
done by considering analytical as well as numerical calculations. The instabilities that will be 
observed will then be treated by using a volume diffusion term that is derived from turbulent 
Reynolds-averaged considerations and implemented into the continuity equation. This term 
is related to an approximate Riemann solver proposed by Ferrari et al. [7j and adapted to 
the present wall boundary conditions. 

The next two sections deal with improvements regarding the imposition of boundary values. 
In Section [5] a new formula for imposing a non-constant driving force based on a volume flux 
will be introduced and compared to a standard formulation. Subsequently, the method of 
imposing Neumann wall boundary conditions will be generalized to Robin boundary condi- 
tions and arbitrary orders of accuracy. A wave equation with Robin boundary conditions 
will be solved numerically in order to demonstrate the capability of the present model. 
The last theoretical contribution presented in this paper will then deal with minor modifica- 
tions to the volume diffusion term as well as the wall boundary conditions to take external 
forces, e.g. gravity, into account. These two modifications are validated by two still water 
simulations. 



Finally, the present methodology is applied to a violent free-surface flow, a schematic dam- 
break over a wedge. Quantitative comparison of wall pressure forces is made with respect 
to other numerical formulations including a Volume-of-Fluid method and it is shown that 
the volume diffusion term avoids using any free-surface correction in the continuity equation 
required in previous work |3]. 

2. SPH Overview 

Smoothed particle hydrodynamics (SPH) is now a well-known Lagrangian numerical 
method. It will be assumed that the reader is familiar with the basics (for theoretical 
grounds refer to Monaghan P] or Violeau [2]). Particles will be denoted by a or 6, and the 
interpolating SPH kernel by w. It is assumed to be radial, i.e. isotropic. The kernel used 
throughout the paper is the S^'^-order Wendland kernel given by 

{ ^(l^lL) (i + 2y\ if0<r//i<2, I,. 

w{r) = l h"'\2h)\h) - / - ) (1) 

[ if 2 < r/h, 

where r is the distance between two particles, m is the dimension, a„i = 7/{Ait) in 2-D and 
h is the smoothing length. The latter in the following is set as /i = 2Ar, where Ar is the 
initial particle spacing. Only weakly compressible SPH [S] in 2-D will be considered in this 
paper. 

Semi-analytical wall boundary conditions for arbitrary boundaries in SPH were introduced 
by Ferrand et al. in [3l [9l [10] . The approach is inspired by the papers from Kulasegaram et 
al. [TT] and Di Monaco et al. [12]. First, the main principles of this method are recalled. 
The key idea is to use an analytical kernel wall correction factor 7 defined by 

7a = / Wafedr^, (2) 

Jo. 

where Wah = ui{\\r^ — rb\\), La being the position of particle a. Here, r^ is used as a continuous 
variable, integrated over the fluid domain Q. If no part of the wall is inside the kernel support 
of a particle a then the value of 7^ is equal to one. 

The boundary conditions are implemented by using the sets of elements described in Table 
[Tj As usual, the fluid is discretized using particles with a spacing initially set to Ar. The 
boundary is split into line segments s of length approximately Ar. These segments are called 
boundary elements and are located between two vertex particles v as shown in Figure [TJ 
Herein, differential operators, both continuous and discrete, will be written in bold {e.g. 
Div). Vectors and matrices will be written as 5 and M respectively. The general discrete 
SPH approximation for a scalar / at position a is given by 

[f]a = -T^VbhWab, (3) 

where H is the volume of particle 6, V the set of all particles and /{, the value of the fleld 
/ associated with particle h. The continuous analog will be denoted with [/]^. Note that in 




Figure 1: The different types of elements. 
Set I Description 



V 
V 

J" 

s 



All particles (V U J^) 

Particles on wall boundary (vertex 

particles) 

Particles inside the fluid domain 

Wall boundary segments 

Table 1: Overview of elements. 



general [f]a 7^ /a, ^-e. the SPH interpolation violates the Kronecker delta property even in 

the continuous case. 

To obtain an approximation of 7, a governing equation is used, given by 



d7a 
At 



Vl ■ Yla, 



where 'tf is the velocity relative to the wall and 



Y7a 



WabUb^b- 



(4) 



(5) 



an 



The normal n^ in the equation above is oriented inwards, a convention used throughout this 
paper. Note that the integral in Eq. (pi) is obtained by a discretization in time and values 
based on the boundary. It is thus distinctly different to the Shepard filter given as 



"a = y^VhWab, 



(6) 



bev 



which originates from a discretization in space based on the fluid particles. 
The other key idea here is to derive rigorously the differential operators without neglecting 
boundary terms coming from the integration by parts as demonstrated below starting from 
the continuous SPH approximation (denoted by the superscript c) of a gradient of a scalar 



field /: 



[Y/]a = - fm)bWa,Ar, (7) 

7a Jn 

1 /■ „_ , 1 



/ fbYb^abdTb H / Yb{fbWab)drh 

Jn la Jn 

/ fbYaWabdTf^ / hWablh^b- 

Jn la Jdn 



la Jn la Jdn 

To obtain the last line the kernel gradient asymmetry (a result of kernel isotropy) was used 
as well as Green's Theorem. From the SPH approximation as given by Eq. ([3]) Ferrand et 
al. proposed the following SPH approximation for the divergence of a vector field 5: 

Div^'^(S) := - — J^m^S,, ■ V^w^b + -$^5,, ■ V7a., (8) 

lapa ^^^ la ^^^ 

where the superscript F stands for "Ferrand" and -B^^ = i?„ — B^ and 

Y.las = / Wablh^b^ (9) 



is the integral of the kernel over an individual segment s, i. e. the contribution of this segment 
to V7a. Hence, the gradient of 7a can be written as 



V7a = 5^Y7a.. (10) 



ses 



Note that the integral in Eq. ^ is known analytically for a given kernel (see Ferrand et al. 
[3] for the case of the Wendland kernel). Following the same paper, the discrete gradient of 
a scalar field / is given by 

Gradr(/) = ^^m, (^ + fj V^w., - f^E (4 + 4) P-V,.,. (11) 

/a ;,gp \Pa Pb / ^"- seS ^^^ ^^^ 

The boundary terms in the above differential operators originate from the fact that the 
differential operator is shifted from the unknown function to the kernel via partial integration. 
In most SPH publications they are neglected. This is not the case in the semi-analytical wall 
boundary conditions where they are transformed to boundary integrals via Green's theorem 
as shown in Eq. ^. 
The continuity equation used in [3J is not the conventional time-dependent form given as 

^ = -PaBW^iv), (12) 

but the equivalent time-independent form 

dilapa) = d f ^rribWab j • (13) 

\bev ) 



In the case of a free-surface flow 7^ is replaced with an heuristic blending function which is a 
necessary correction as otherwise particle repulsion can be observed as demonstrated later. 
The Navier-Stokes equation is given as 

^ = -lGrad^'^(p) + Lap^ (/x, t;) + g_. (14) 

Cl'' Pa 



where the SPH Laplacian is defined below by Eq. (18). The system of equations is closed 
with the Tait equation of state [13J , given as 



Pocl 



Po 



(15) 



where po is the reference density, cq the numerical speed of sound and ( = 7 for water. 
Continuing to follow [3] , in order to impose Neumann boundary conditions at the wall for a 
scalar / the values from the fluid can be extrapolated to a vertex particle v via 



fv = — y^Vb fb w^b, (16) 

where 



"^^b^T 



a^ = ^Vb Wbv (17) 

The values for the boundary segments s are given as an arithmetic mean of the associated 

vertex particles. The above represents a first-order approximation of the true Neumann 

condition. 

Finally, some aspects of the viscosity treatment as proposed in [3J will be summarized. The 

first requirement is a Laplacian operator for which the same approach as Morris et al. |14j 

is used but including the boundary terms. Ferrand et al. [3J approximate the Laplacian as 

V-(/V®5) ^ Lap^'^(/,5) = P^Y.mJ-^^±^^r_,,.V,Wa,--Y,fs{^B)s-Y-ias, (18) 

^<^ feg-p PaPb ^ah 7a ^^^ 

where Vah = WLabW- ^^ ^^e case of viscosity B = v,f = p and so /is(V (g) t^)^ ■ z^^ = r^, i. e. the 
wall shear stress vector at segment s. The latter is given by 

dv 

PQ^=Phr\kr, (19) 

where, by definition, v^ is the friction velocity. In the case of laminar flow this is approxi- 
mated by 

11^.11^. = ^, (20) 

where 2; is a short distance from the wall and v = p/p. This equation only holds for z '^ L, 
where L is a characteristic length-scale of the flow geometry. It can be extended to turbulent 
flows using wall functions, see [3]. 
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3. On the skew-adjoint property including boundaries 

3.1. Theoretical investigation 

Now, new developments will be presented in order to better understand and improve 
Ferrand et al.^s [3] model. 

In this section the focus will be on the property of skew-adjointness of the two arbitrary 
(discrete or continuous) operators Grad and Div. These two operators are skew-adjoint if 
and only if 

< Grad(/), 5 >=-</, Div(5) >, (21) 

where () are the respective L^ scalar products. Before starting the actual investigation the 
importance of this property shall be highlighted. Consider a system of particles representing 
a fluid without external influence nor dissipative forces, then its energy is given by 






E = Ekin + Eint = 2_^ma [ -\\Va\\ + eint,a{pa) ] , (22) 

where E^in and Eint are the total kinetic and internal energy respectively, while eint,a is the 
specific internal energy of particle a. The time derivative of the kinetic energy is given by 

—^ = X]rnat;„ ■ -^ = -J2^aVa " — Grada(p), (23) 



where the Euler momentum equation was used, i. e. Eq. ( 14 ) without the Lap and g terms. 
The time derivative of the internal energy can in turn be written as 

where the last equality follows from thermodynamic principles that relate the internal energy 
per unit mass to pressure and density via p = p^ dcint/dp. Looking at the time derivative of 
the total energy one thus obtains with the help of the continuity equation (Eq. (12)): 

^ = -Yy-"^- ■ Grad,(p) - 5^KPaDiv.(t;). (25) 

Written in notation with discrete scalar products this yields 

^ = - (Grad,(p),t;J - (p„Div„(t;)) . (26) 

This shows that the energy is exactly conserved if the two discrete differential operators 
Grad and Div are skew-adjoint (Eq. ([2T|), i.e. if the right-hand-side of (26) is equal to 
zero. This property is natural, since the same occurs with ordinary differential operators 
when changing the discrete sums with integrals (see e.g. Violeau |2j): 

Skew-adjoint operator definition: (V/, 5) + (/, V ■ 5) = — / f{r)B{r)-n{r)dr. (27) 

^ V ' Jan 

= :SA 



li f = p and B_ = v, the pressure and velocity, respectively, which solve the Navier-Stokes 
equations, then the right-hand side is equal to zero. This is due to p = at a free-surface 
and w ■ n = at a solid wall. It would thus be advantageous if the SPH operators adhere to 
this property. 

The investigation into skew-adjoint operators will commence by first considering basic con- 
tinuous SPH operators without boundary terms, i.e. using integrals instead of discrete sums 
and assuming that no boundaries are present. The operators under investigation are given 
by 

Grad^'^(/) = [ hYaWabdr, (28) 

Jn 

Div^'^(5) = f B,-YaWabCk,, (29) 

JQ 

where the superscript b stands for "basic" in the sense that it is the SPH operator which 
can be derived directly from the interpolation without the addition of any other terms. The 
superscript c stands for "continuous" similar to the SPH interpolation. The left-hand side 



of the skew-adjoint operator definition (Eq. (27)) is then given by 



SA= f f ih V^w^t -B^ + faR,- V^wab) dr,dr,. (30) 

Jn Jq 

Due to integral additivity it is possible to swap the dummy labels a, b in the second term 
which yields 



SA= ih B^ ■ YaWab + hB^- Y.kWab) dr.dr,. (31) 

Jn Jn 

The kernel isotropy result in the asymmetry of its gradient, i.e. V_jVab = — V^'U^afe, shows 

finally that SA = 0, i.e. that Grad ''^ and Div '^ are skew-adjoint. Note that this property 

also holds true in the discrete case, where the integrals are replaced by sums, from the same 

properties. 

Now the following symmetrised operators shall be defined: 

is,c,kfr\ _ f Pa h + Pb fa. 



Jn PaPh 

Divr'^(5) = 4 / p^p^(5, - 5J ■ V,w.,dr„ (33) 

Pa Jn 

where the superscript s stands for "standard" and fc is a power used to discuss a wide range 
of operators using a sole notation. Its effect will be discussed in the following. Note that. 



apart from the lack of boundary terms, the symmetrized gradient (32) differs from the dis- 



crete form (11) by the presence of density terms. Eq. (11) can be recovered by setting k = 1. 



The above operators in Equations (32) and (33) are skew-adjoint if the newly added terms 



have opposing signs. The proof, which will be omitted, can be found e.g. in [2J. As in the 



previous case it holds true in the discrete case also. 

When removing the assumption of no boundaries the calculations are no longer as straight- 
forward. Following a similar procedure to Equation ([T]), in the vicinity of a boundary the 
operators of interest are given in continuous form by 



Gradr''(/) 



BWr''{B) 



1 f PfAlPf^ v...d.. 

la Jn PaPb 

la Jan PaPb 

1 



(34) 



yaPa Jn 

- ^ PaPtiRb - B^) ■ nf^Wabdr^ 

laPa Jan 



(35) 



where the superscript 7 indicates the renormalization. Apart from the additional density 
terms in the denominator of the Grad operator, these formulae are the continuous analog 
to Eqs. (Sl) and (11) for k = 1. In Appendix A using the gradient and divergence given by 
Eqs. (34) and (35) the left-hand side of the skew-adjoint operator definition, Eq. (27), can 



be shown to reduce to 



SA 



(h^O) 



falL ■ Tka^a 



(36) 



an 



This is equivalent to Eq. (27) showing the skew-adjoint property for SPH continuous oper- 



ators with boundary terms in case of the correct imposition of the boundary conditions in 
the limit /i — )■ 0. Note that it is essential for this result that the operators are renormalized 
with 7. Furthermore, the process of taking the limit is only required due to the violation of 
the Kronecker delta property by the SPH interpolation. 

The derivation shown in [Appendix A could have been made significantly shorter, if we were 
only interested in continuous operators. However, the reason for explicitly going through this 
derivation was to show the steps necessary to prove the same in the case of discrete operators. 
It would thus be necessary to have a discrete version of the Stokes' theorem. Additionally, 
the Kronecker delta property is required, but that is violated even by the continuous SPH 
interpolation. The difficulty for a discrete Stokes' theorem is due to the fact that V7 is 
calculated analytically along the wall (see Section [2]), whereas the volumetric integral over 
V_w is approximated via a discrete sum. 

As a consequence, discrete SPH operators with boundary terms as presented here are not ex- 
actly skew-adjoint contrary to continuous SPH operators. Inside the fluid, however, (7 = 1, 
no boundary terms) the skew-adjoint property is fulfilled in the discrete case as it is equiva- 
lent to the standard differential operators (Div'*''^''^, Grad*''^''^) given by Eqs. (32) and (33). 



Finally, the discretized forms of Eqs. (34) and (35) are given by 



Grad-^/) = -E^ /" ^'t'" ^" V.^^^- -E """ ^'t'' ^" V7.. (37) 

laPa f,gp laPa ^^g 

Note that these operators are a generalization from the ones derived by Ferrand et al. given 



by Eqs. dSl) and (11) (where k = 1). In the following fc = will be used unless otherwise 
noted. According to the authors experience this choice has at best a marginal influence on 
the results. 

3.2. A numerical experiment 

As the above theoretical investigation indicates the SPH discrete operators for the semi- 
analytical boundary conditions are not skew-adjoint, and thus do not conserve total energy. 
Hence, a flow between two inflnite plates (at 2; = and z = h) will be used to address this 



issue from a numerical perspective. Equations (12) and n\.A\ are solved with the pressure 



calculated via the equation of state ( 15 ) and the particles are moved according to dr^/dt = v^ 



The viscosity is set to zero and the following velocity proflle is used as initial condition 

Vx = 0,t;^ = -^sin(47rz), 

where Cq is the numerical speed of sound and z G [0, 1]. The quantity of interest is the energy 
E which is deflned by 

t 
E{t) = ^(<Grad^(p),t;> + <p,Div^(t;) >)" , (39) 

n=l 

where the subscript n indicates the current time-step. The value of E is equal to the time 



integral of the left hand side of Eq. (27) and it is zero if the operators under investigation 
are skew- adjoint. 

In Figure 0the non-dimensionalized value E^ = E /{pQC^h?) can be seen plotted over time for 
the present boundary condition as well as the ghost particle approach [15]. Two important 
features of this plot shall be highlighted. The flrst is that the value E is never zero, indicating 
that the operators are indeed not analytically skew-adjoint in case where boundaries are 
present. Compared to the ghost particle formulation the energy introduced to the flow is 
lower with the present boundary conditions. 

It shall be noted that this simulation was also run with the formulation by Kulasegaram [TT] 
which resembles the present approach but, in contrast to the present method, is analytically 
skew-adjoint. However, as Ferrand [3J already pointed out, the operators by Kulasegaram 
lack a term in the gradient and thus do not represent the physics accurately. The simulation 
was also conducted with the formulation by Monaghan and Kajtar [5j but the results were 
completely non-physical. 
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Figure 2: Non-dimensional energy budget over time i"*" = tco/h for ghost particles and present approach 
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Figure 3: Values of v+ at different times {t+ = 0.0: Fig. 3(a) t+ = 9.0: Fig. 3(b) ) 



The second observation is that the flow starts exhibiting non-physical fluctuations after some 
time. In Figure [sjthe flow can be seen at t"*" = tco/h = and t+ = 9.0 (corresponding to 
the times indicated by black circles in Figure [2]) where the non-dimensionalized values of v is 
plotted at each particle. It was already remarked in the previous section that the deficiency 
has to originate from the boundaries, as inside the fluid the operators are skew-adjoint. 
This result can be connected to the numerical experiments of Basa et al. [16j who simulate 
the laminar Poiseuille flow with viscosity. They also observe that the break-up that occurs, 
originates from the boundaries as particles which are close to the wall experience different 
viscous forces due to numerical fluctuations. It should be noted that the behaviour of the 
viscous system is quite similar to that of a laminar flow that becomes turbulent. Starting 
from this observation the continuity equation will be investigated in the next section to 
extract the mean flow field in order to suppress this numerical "turbulence" . 
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4. A volume diffusion term for numerical turbulence 

4.I. Basic idea 

In the context of turbulent flows the continuity equation in the Reynolds averaged context 
is given by 

J- = -pY-v-Y-i7v^), (40) 

dt 
where the primes refer to turbulent fluctuations and the bars the Reynolds averaging (see 
e.g. Pope |l1^). Moreover, the bar on the Lagrangian derivative means that the advection 
velocity is in principle the averaged one. In laminar flows, as considered here, the fluctuating 



quantities are supposed to be zero. Still Eq. (40) can be applied to numerical fluctuations 



in an attempt to stabilize them. The gradient- diffusion hypothesis states that 



p'v; = -KYp, (41) 

where K is the turbulent diffusivity. Inserting this into the averaged continuity equation 
yields 

S- = -pV-v + Y-{KYp)- (42) 

dt 

As no fluctuating quantities remain we will drop the bars in the following. Expressing Eq. 

flowing is obtained 

p„Div^(t;) + Lap^(ir,p), (43) 



(42) in terms of SPH operators the following is obtained 

dpa 
dt 

where Lap''' is an SPH discrete operator, here applied to p with a diffusion coefficient K. If 



the model by Morris [lAl is used for the discretisation of the Lap'*' operator as in Eq. (18), 
without boundary terms, then the full discretisation reads 

^ = Pa J]H U, + {Ka + Kt)pab^) " Ya^ab, (44) 

where pab = Pa — Pb- It is common to write the diffusivity term as -ft' = vt/<^t-, where vt 
is the turbulent viscosity and aT is the turbulent Prandtl number. Continuing the analogy 
with physically-based turbulence, one may use the mixing length model to estimate vt (see 
e.g. Pope [Ej), leading to 

K^-Llj, (45) 

where U and L are characteristic velocity and length scales, respectively, and L^ is the 
mixing length. Deflning the numerical Mach number as Ma = U/cq, where again Cq is the 

numerical speed of sound, yields 

Ma 2 Co , . 

K ~ -—L^-. (46) 



12 



Typically, ctt ~ 1, I^^ 



L/10 and in weakly compressible SPH Ma = 1/10 [T], which yields 
ir~A^. (47) 



Ar 103 
Depending on the resolution, K can thus be given as 

Co Ar 



K 



V 



(48) 



where rj ^ lO^Ar/L typically has the range of values of (9(10) — (9(100) depending on the 

ratio of Ar/L used to resolve the length scale L. 

Ferrari et al. [7j proposed a diffusion term which is remarkably similar to the one above (Eqs. 



(44) and (48)). It is based on the theory of Riemann solvers which were first introduced to 
Ferrari et al. use an approximate Riemann solver to obtain the following 



SPH by Vila 
continuity equation 



where 



dpa 
dt 



E^^ 



feeP 



Lab 



h I PbVab + Cab—^Pab 
Tab 



Cab = max(Ca,Cfe), 



YaWab, 



and 




(49) 



(50) 



(51) 



with po being the reference density and ( the exponent of the equation of state (15). Com- 



paring Eqs. (44) and (49) it can be deduced that our model is equivalent to that of Ferrari 
et al. [7] if 

rabCab = Ka + K^. (52) 



Now a closer examination of Ca, given by (51), Eq. (15) shows that 



Co 



c-i 

2 




(53) 



so Ca is nothing else than the speed of sound of particle a. Thus K has the dimension m'^/s 
as expected for the turbulent diffusivity term. 



Comparing Eqs. (52) and (48) shows that the correction as proposed by Ferrari et al. [7j 
applies a significantly higher viscosity term than what can be derived from the simple mixing 
length model. This fact can be of importance when it comes to the use of this correction 
in the context of turbulent flows, where SPH is used increasingly {e.g. Ting et al. [19]). 
The turbulent viscosity introduced due to the volume diffusion term by Ferrari et al. [7] 
can dissipate more than just the numerical noise and thus have an influence on the energy 
spectrum of the flow. 

Finally, note that in the Reynolds averaged context the Navier-Stokes equations would also 
need to be averaged. This was neglected in the above as the aim was to find a different 
interpretation for the volume diffusion term by Ferrari et al. [7] which only acts on the 
continuity equation. 
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4-2. Semi- analytical wall boundary framework 

After this interpretation the question arises on how this additional numerical diffusion 
term can be included in the wall boundary formulation as described above. To do so, the 
flux of the quantity KVp has to be investigated in normal direction of the wall. As the 
boundary condition on the density implies dp/dn = for flows without external forces, this 
flux is zero as well. Thus, using the Laplacian of Ferrand [3], |9] and the divergence given by 



Eq. (38) the continuity equation with volume diffusion term reads 



-E^^ U, + -Pa,"=^) ■ YaW^, - ^-Y.V^s ■ Y.las. (54) 

I , ,-73 V Pa ^ab / 7a „, c 



'^ _ PaS^^r ( .. , '^«'' „ ^ab\ w... P' 

dt ~ % 



"■ bev ^ '""■ '""^ '"se5 



As can be seen from the above, the volumetric term contains also vertex particles. In turn 
however, their density is determined from the boundary condition. As this volume diffusion 
term can be seen as a way of transferring volume from one particle to the next it is important 
that this term is anti-symmetric, i. e. the volume taken from a particle a due to the influence 
of particle b needs to be added to particle 6 as a result of the influence of particle a. This 
principle is thus violated if vertex particles are taken into account in the volumetric sum. In 
2-D simulations this change has negligible impact on the flow. 

5. Imposing a volume flux in periodic viscous flows 

Related to the issue of solid boundaries are open boundary conditions. In the following 
the focus will lie on periodic boundaries which itself do not pose an issue with SPH. However, 
in order to drive such a periodic flow with a fixed volume flux and thus variable force only 
one formula exists, which shows relatively large errors as will be shown in the following. 
Imposing movement in the standard Poiseuille flow is normally achieved with a fixed driving 
force. This is however only possible as the value of the latter was known a priori, since the 
expected velocity profile is known. It is more common to drive a flow with a certain volume 
flux Q. In the following it will be described how to impose a variable driving force based on 
the expected volume flux. 
The volume flux Q is defined as 



Q = / v-dA= / V- ndA, (55) 

Ja J a 

where A is the cross-section area of the flow, n its normal and v the velocity. To obtain this 
value in the SPH framework an average over all particles in a slice of the domain will be 
taken, i.e. 

QsPH = -^^ 5Z ^^'^b ' ^^' ^^^^ 

where Ar^ is the width of the slice and the set Va contains the particles in the slice. Herein, 
Ava is twice the particle spacing Ar. This can also be written as 

QsPH = Av, (57) 
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where v is the cross-averaged velocity. 

In the book by Violeau [2] this average velocity is used to calculate the force via 



.e.t,n V~2v--'+V-~^ 



2At 



{5t 



where n is the n-th time-step, F is the force in direction of the normal and v = Q/A the 
desired velocity. This formula originates from the finite volume community [2nj- As can be 
seen in Figure |4] this value does not converge to the analytical (expected) one. The reason 
for this deficiency comes from the fact that internal forces are not considered in the above 
formulation and so the external force is always underestimated. The equilibrium that should 
be reached is defined by 

pext ^ pint ^ g^ ^gg^ 

where F^^* and F™* are the external and internal force respectively (the latter includes 
pressure and viscous forces). If the system is not in equilibrium the following holds 



nex4,n 



+ F 



int.n 



fj y'' 1 1 



Ar- 



(60) 



Ideally the velocity reached at time n is equal to the desired velocity v = Q/A. So rewriting 
the above yields 



7:ext,n 



V — V 



n-1 



;^int,n 



At" 



(61) 
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Figure 4: Comparison of the error in the bulk velocity of the Poiseuille flow. 
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Clearly _F*"*'" is not available but it can be assumed to vary only little between two consec- 
utive time-steps, i.e. _p*"-*'"- ^ pmt,n-i^ The latter value can then be calculated using Eq. 



(60) to give 

^ At" At"-i ■ ^ ' 

Finally rearranging the above yields the following new formula for calculating the external 
force: 

7; _ 9?T"-1 -I- ^"-2 
Tpext,n __ ^ '_^ |_ Tpext,n-1 (fyK) 

where it is assumed that the time-step At is constant. The first term in the right-hand 



side of Eq. ( [63| is twice that of Eq. (58). Thus (58) amounts to considering that F^^* is 
independent of time, which is a crude approximation. 

In Figure |4] the two different means of imposing a driving force based on a volume flux are 
compared in the Poiseuille flow case. To do so the relative error in the bulk velocity is 
plotted over time. In Figure |4] it can be seen that with this new formulation the external 
force converges much more closely to the theoretical value. In the case of the Poiseuille flow 
it is not necessary to impose a volume flux but instead the analytical force can be used. 
However, in cases where the analytical value of the internal (viscous) force is not known a 
priori the above formulation provides the means to impose a volume flux which converges 
to the desired value. 

As it can be seen the present approach reduced the error by about flve orders of magnitude. 
The error obtained with the original formulation is close to 1% which is not negligible. 

6. Generalization of wall boundary conditions 

After this analysis of the unifled semi-analytical wall boundary conditions and the rein- 
terpretation of the volume diffusion term the focus in the following two sections will shift 
towards novel developments which expand the boundary model. 

6.1. Theory 

As mentioned in Section [21 to satisfy von Neumann boundary conditions Ferrand et al. 



[3] showed a first-order approximation approach (see Eqs. (16) and (17)). In the following 
this approach shall be generalized to arbitrary orders of accuracy and to Robin boundary 
conditions, which include Neumann boundary conditions as particular case. Ryan et al. [21j 
implemented Robin boundary conditions for SPH with the use of an additional source term 
in the governing equation. In this section a different approach will be shown that enables 
Robin boundary conditions to be imposed directly. 
Such a boundary condition is given by 






= /is, (64) 
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for an arbitrary scalar field / with given /ii,yU2 7^ 0,/i3. To impose Neumann boundary 
conditions (/ii,/i2) take the values (0, 1). The above can be rewritten as 



dn /i2 



(65) 



The key idea is to use a weighted linear least squares approximation of the desired field in 
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Figure 5: Illustration of the arbitrary order Robin boundary conditions for SPH. Blue particles denote fluid 
particles and particle v in orange is a vertex particle under consideration. 

the direction of the normal n^ of a vertex particle f , which is defined as the average over the 
two adjacent segment normals. This implies that the problem is considered locally as one 
dimensional projected onto the normal n^ of the vertex particle v. The procedure presented 
in detail in the following is illustrated in Fig. [5j 
Let m be the desired order of approximation and define a polynomial A as 



A(x) 



PiX H x + l3i, 



i=2 



1^2 



(66) 



X being the normal distance to the wall. 

We set (3i = flg^, i^ (3i is the value of / prescribed at the wall. The polynomial can also 

be written as 

— X + Pi 1 ; 

/^2 V /^2 



A(x) 



^Ax' + —X + I3i { 1 - —X 



4 = 2 



(67) 



This represents an approximation of / in the direction of the normal n^ with /^ = A(0). For 
a fiuid particle a E J^ and Xa = Hav ' Tkv the value of A^ is thus given as A(r^^ ■ n^) when no 
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external forces are present. Let X G M(|J-'|,m) be a matrix defined as 

X =[ 1-S(^--^^) if^' = l' (68) 

I J-"! is the number of fluid particles and defining y E M'-^I as 

ya = fa -{La^-n^)- (69) 

The goal is to minimize fa — Xa- Ifence, we are looking for a solution of the over constrained 
system of linear equations 

K-P = y, (70) 

where /3 is the vector of components /3j and the dot represents a single contraction, in this 
case the matrix multiplication. A common approach to this type of over-constrained problem 
is to use the least squares method which can be solved via 

K^ ■K-/^ = K^ -y- (71) 

In order to put more weight on particles closer to the vertex particle, the SPH kernel is used 
to construct a weighted least squares interpolation. To include this into the above let Sab be 
the Kronecker delta and A G MdJ-"!), a square matrix, be defined with elements given by 

=ab ^ ^abVaWya, (72) 

(there is no summation over indices here). Then the weighted least squares interpolation 
can be found by solving 

K^ ■A-K-l3_ = K^ -A-y- (73) 

Due to the fact that the problem is one dimensional the matrix on the left-hand side is of 
size m X m and can be inverted easily. Of interest is A(0) which is /3i. Thus finally 

/, = A(0) = /3i = ((^^ ■ 4 ■ K)-' ■K'-t M)i' (74) 

where the subscript 1 on the right hand side refers to the first component of the vector in 
brackets. 



If m = 1 the above reduces to the equation shown in Eq. (16) thus showing that this 
indeed is a generalization to arbitrary order. Note that the matrix that needs to be inverted 
{X ■ A • X_) is contained in M(m) and thus to obtain second-order boundary conditions the 
inversion of a 2 x 2 matrix is required in two dimensions. 

Finally, two things should be noted. Firstly, the matrix is non-degenerate if at least m fluid 
particles are in the neighbourhood of v and that they have different values of (r„^ ■ n^). 
Secondly, the formulation as presented above is independent of the method used to describe 
the walls, so it could also be used e.g. for the SPH ghost particle approach [6] . 



18 



Analytical solution 
Fist order BCs 
Second order BCs 
^Neumann boundaries 




Analytical solution 
Fist order BCs 
Second order BCs 
[""-., Neumann boundaries 




(b) t = 0.22 



0.5 



-0.5 



Analytical solution 

Fist order BCs 

""-■--Second order BCs 

Neurfiann boundaries 




(c) t = 0.7 



0.5 



-0.5 



Analytical solution 

Fist order BCs 

Second order BCs 

--., Neumann boundaries 




0.2 0.4 0.6 O.i 

X 



(d) t = 0.81 



Figure 6: Wave equation (75 1 with Robin boundary conditions (76) at four different time-steps. 



6.2. The wave equation as an example 

In this section it will be shown that the above formulation indeed works, i.e. that Robin 
boundary conditions can be enforced. For this purpose, a wave equation is considered in 1 
dimension: 



with boundary conditions 



du 
dn 



df^ dx"^ ' 



({0,l},t)=M({0,l},t), 



and initial conditions 



u{x, 0) = and 



du 



X, 0) = sin(27rx). 



(76) 



(77) 



The domain [0, 1] is discretized with an initial particle distance Ar = 0.01 and the particles 
are fixed. The second-order spatial differential operator is again discretized by using the 
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Laplacian as given in Eq. (18). The temporal derivative is discretized by a second-order 
finite difference scheme. The analytical solution [22] of this problem as function of the 
dimensionless variables x and t is given by 



Ua 



{x, t) = y^Oii sin(Kjt) [sm{K,ix) + Ki cos(fi;jx)] 



(7J 



where Hi is given as the i-th root of 



iaii{ni 



ZKi 



KJ 



(79) 



and the ctj are uniquely determined by the initial condition. 

First, a qualitative view on the solution is given in Fig. [6] at times t = {0.1, 0.22, 0.7, 0.82}. 
It can be seen that the main features of the field u are reproduced by the SPH solution and 
that, due to the new boundary condition formulation, the Robin boundary conditions are 
correctly imposed as we observed a distinct difference in the solution when using Neumann 
boundary conditions, i.e. /ii = 0. In order to illustrate this the analytical solution for the 
Neumann boundary case was superimposed in Fig. |6J The snapshots in Figs. 6(a) and 



pi The snapshots in Figs. 
6(d)| are plotted at instants where the error on the boundary is maximal. Compared to 
these snapshots, the error between the SPH solution and the analytical solution is smaller 



throughout the domain at intermediate times (Figs. 6(b) and 6(c)) 
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Figure 7: Wave equation with Robin boundary conditions: Relative error of the SPH solution over time at 
the boundary. 



There is a notable difference in the SPH solution depending on whether the boundary con- 



ditions are of first or second order, i.e. whether ?ti = 1 or m = 2 in Eq. (66). To quantify 
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this more precisely consider Fig. [7] where the error at the boundary is plotted over time. It 
can be observed that the error is sensitive to the order of the boundary condition and that it 
can be reduced by up to 30% by choosing a second-order approximation. The vertical bars 
in Fig. [7] show the instants that were shown in Fig. |6| 

7. Still water with a free-surface 

In the following, two of the above presented improvements will be reviewed in the context 
of still water free-surface evolution. The first section below will focus on the volume diffusion 
term and the second one on the arbitrary order boundary conditions. 

7.1. Modification of the volume diffusion term 

Using the volume diffusion term presented in Section |4] in a simulation of still water shows 
that the term as shown above will not have a zero contribution towards the density equation. 



This is caused by the fact that the boundary terms in Eq. (54) do not vanish when gravity 



is present. However, there are no segments on the free-surface and the problem can thus not 

be resolved by adding a boundary term. 

In order to compensate for this deficiency the correction can be modified so that it reads 



dpa 
dt 



Pa 

la 



Eh 

bev 



I Cab Lab 
ab "T Qab 

Pa Tab 



■ YaWab - — y'l^a. " Yla 



(80) 



in place of Eq. (54) where the density difference pab was replaced by the modified density 



difference Qat which is given by 



Qab — Pa — Pb 



gpof 

-^{Zb 
'-0 



Za)- 



(81) 



where g is gravitational constant and z the vertical elevation. This is a linear approximation 
of the difference between the non-hydrostatic densities of two particles. The external force 
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Figure 8: Position of free-surface in infinite open channel at rest. 
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in the formula above is gravity, but adaptation to other forces is straightforward. 



In Section |4.2| it was remarked that the volume diffusion term should be anti-symmetric with 
respect to particles a and h in order to avoid that the global volume is changed. Clearly Eq. 
(81) obeys this principle. A similar correction was simultaneously proposed by Antuono et 
23|. 



al ^ ^ 

To illustrate the difference between using the traditional pab = Pa—Pb and Qab consider Fig. [8} 
which shows the non-dimensionalized position of the free-surface of an open channel-flow at 
rest where the initial density is set to po- In the plot the time was renormalized by hgwi/co, 
where hswi = 1 is the still water-level and cq = 10 the numerical speed of sound. It can 
be seen that without the above modification the free-surface detaches due to a transfer of 
volume from the denser lower part to the upper part of the fluid. The modification proposed 
above clearly avoids this issue keeping the free-surface elevation almost constant as expected. 
The decrease is due to the initial condition and the weak compressibility of the fluid. 

7.2. Modification of the boundary interpolation 

In this section the generalized boundary conditions presented in Section [6] will be investi- 
gated in the presence of gravity. This means that the function / will be equal to the pressure 
for which the classical Neumann boundary condition {dp/dn = pg ■ n) will be applied. 
The approach presented in Section [6] will produce unsatisfactory results at the intersection 
of free-surface and a wall as tangential variations are not neglected. The constraint above 



(Eq. (69)) reads 



A. 



Va = Pa- 



(82) 



In order to neglect tangential variations for external forces such as gravity, the proper con- 
straint is given by 



K^ya=Pa- Pa 
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Figure 9: Still water in a closed tank (left: without correction (82), right: with correction (83)) 
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In Fig. [olthe difference between Eqs. (82) and (83) is shown. The resolution is chosen to be 



relatively low in order to highlight the impact of the proposed correction. The picture on the 
left hand side shows particle movement which is an order of magnitude larger than the one on 
the right hand side which demonstrates the corrected approximation. Similar to the velocity 
field, the pressure prediction is improved as well by lowering the magnitude of pressure waves 
originating from this corner. To explain the formula presented above consider the setup in 
Fig. |9] with a perfect hydrostatic pressure distribution. Now we look at a vertex particle on 
a vertical wall which is located next to the free-surface. When constructing the polynomial 
A (Eq. ([66))) as given in Section [6] the fluid particles considered for the approximation all 
have a pressure greater or equal to zero. This causes the pressure of the vertex particle to 
be greater than zero, although its theoretical value is zero. This in turn causes a repulsive 
force that can be seen on the left hand side of Fig. |9| If, on the other hand, the hydrostatic 



part is subtracted from the fluid particles as in Eq. (83) then all fluid particles used for 



the approximation of A will have zero pressure and thus the vertex particle has the correct 
pressure. 

8. Dam-break with wedge 
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Figure 10: Schematic dam-break on a wedge: Comparison of forces on left wedge wall. 

After having analyzed and extended the present wall boundary conditions a final simulation 
shall be performed. This uses most of the theoretical results presented above in a more 
complex free-surface flow. 
A schematic dam-break over a wedge will be simulated with the same geometry as used in 
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Figure 11: Comparison between VOF and SPH. 
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Figure 12: Steady flow on a wedge: Without (Fig. 12(a)[ ) and with (Fig. 12(b) ) the volume diffusion term 
(84) to avoid particle repulsion. 



[3]. The initial volume of water is 1 in high and 0.5 m wide and the dynamic viscosity is 
set to 1/ = O.Olm^/s resulting in a Reynolds number of approximately 140 based on the 
maximum velocity. The box has a length of 2.2 m where the wedge begins after 0.85 m with 
a side-length of 0.25 m. The force will be calculated as the integral of the pressure along the 
left wedge wall. A Volume-of-Fluid (VOF) simulation is taken as reference solution (Open- 
Foam [21] )• It should be noted that the latter is a multiphase simulation and thus some 
discrepancies are to be expected when compared to the SPH single-phase simulation as il- 



lustrated in Fig. 11, As shown, the traditional boundary conditions using fictitious particles 
[25] or the Lennard- Jones potential [1] fail in predicting the force. Comparing the approach 
by Ferrand et al. [3] with the present one, it can be seen that the volume diffusion term 
successfully reduces the numerical noise, while still showing closer agreement with VOF in 
Fig. 
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Finally, Fig. 12 (a) [ shows the steady state solution of the dam break with the time-independent 
continuity equation (Eq. ([I3])) without the heuristic free-surface correction. Particles that 
were initially on the free surface have retained their larger volumes producing the strange 



bubbles and unphysical pressures. Following the developments in Section 4.1 the integrated- 



in-time volume diffusion term was added, so that the continuity equation now reads 



7>: 
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where the superscripts refer to the time iteration. Note that again the gravitationally cor- 
rected volume diffusion term (Eq. ( [8l| ) is used through gab- The continuity equation above 
is equivalent to using Div''''*^ (Eq. (38)) with k = 1 so for consistency Grad''''^ (Eq. (37)) 



has to be used for the discretisation of the pressure gradient. 



From Fig. 12(b) it can thus be concluded that the heuristic free-surface correction is no 
longer required (contrary to [3|) as the volume diffusion term successfully redistributes the 
higher volume of particles initially on the free surface once they become entrained in the 
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fluid body. 

9. Conclusion 

The semi-analytical wall boundary conditions for SPH introduced by Ferrand et al. [3J 
were recalled in Section [2j The skew-adjoint property of discrete operators was investigated 
theoretically and numerically, showing that it does not hold in the discrete case leading to 
some errors in the conservation of energy. As shown by Morinishi et al. |26j for non-uniform 
grids exact conservation is not required if the error is small. A detailed error analysis of the 
conservation of different boundary conditions would be of interest for further research. An- 
other general issue with SPH is instability due to the existence of zero-energy modes which 
manifest themselves in numerical noise. This noise can be explained by comparison with 
a Reynolds averaged continuity equation which is equivalent to the approximate Riemann 
solver by Ferrari et al. [7]. This interpretation justifies the addition of an additional con- 
stant which was shown to depend on the relative resolution. As the volume diffusion term 
introduces an artificial viscosity this constant prevents excessive damping which would be 
problematic e.g. in the simulation of turbulent flows. 

Subsequently a novel formulation to prescribe a variable driving force for flows with imposed 
volume flux was presented which improves the predicted flow rate by 5 orders of magnitude. 
Additionally, the Neumann boundary conditions by Ferrand et al. [3J were generalized to 
arbitrary orders of interpolation and Robin type boundary conditions. The formulation was 
shown to correctly impose Robin boundary conditions and so extending the possible appli- 
cation areas. Finally, two modiflcations to the boundary conditions as well as the volume 
diffusion term were presented in order to correctly deal with free-surface flows, reducing 
unphysical velocities at the surface by at least an order of magnitude for still water. 
In Section |8] a dam-break over a wedge was simulated showing the capabilities of the present 
improved model compared to a well-known VOF code. This simulation was also used to 
highlight the fact that the volume diffusion term can correct the free surface when using the 
time- independent continuity equation, as proposed by Ferrand et al. |3]. 
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Appendix A. Derivation of skew-adjointness 



In the following the left hand side of Eq. (27) will be analyzed when replacing the nabla 



operator with the gradient and divergence given by Eqs. (34) and (35). Initially the part 



containing the volume integrals of SA Eq. (27) will be investigated. 

>f A + Pl'fa 



O/l,, 




'n Jn la 
1 



PaPl 



'-B^ + fa^-^iB, 



pf 



B„ 



■ YaWab^b^ 



(A.l) 



o 

, , ~lfbBa + —nfaBf^ 

In Jn la LPfe Pa 
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Furthermore, due to the additivity of the integral we can spht it up 



bAr, 




faB, ■ Ya^abdr^dr, (A.2) 
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Jn Jq la Pb 

In the last line the kernel gradient asymmetry V^^afc = — V^'U^ab was used in both terms. 
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The boundary part of Eq. (27) can be reformulated to 
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Combining the volumic SA^ and the boundary term SAb yields 

SA = SA^ + SAb 
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where in the last step a reverse integration by parts was used. The terms in square brackets 
represent SPH approximations which in the limit of /i — )■ converge to 
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where in the second last line again a reverse integration by parts was used and the final line 
follows from Stokes' theorem. 



29 



